#!/usr/bin/perl
#Programmer: Rory Carmichael
#Purpose: Find an appropriate rarity cutoff for a given discretization of a sample and a discretization of the 200 rrts of that sample
#Usage: calc-cutoff.pl <sample-discretization> <rrt-discretization>

if ($#ARGV != 1) {
	die("Usage: calc-cutoff.pl <sample-discretization> <rrt-discretization>\n");
}
my $sampdisc = $ARGV[0];
my $rrtdisc = $ARGV[1];
my @samphash = ();
my @rrthash = ((\@samphash) x 200);
my @rrtmeans = ( (0) x 200 );
my @rrtstdev = ( (0) x 200 );
open(IF, $sampdisc);
foreach my $line (<IF>) {
	chomp $line;
	my @splitline = split(/\t/, $line);
	$samphash[$splitline[0]] = $splitline[1];
}
close(IF);
open(IF, $rrtdisc);
foreach my $line (<IF>) {
	chomp $line;
	if ($line !~ /^=/) {
		my @splitline = split(/\t/, $line);
		push (@{$rtthash[$splitline[0]]}, $splitline[1]);
	}
}
close(IF);
for(my $i=0;$i<200;$i++) {
	foreach my $val (@{$rtthash[$i]}) {
		$rrtmeans[$i] += $val;
	}
	$rrtmeans[$i] /= scalar @{$rtthash[$i]};
}
for(my $i=0;$i<200;$i++) {
	foreach my $val (@{$rtthash[$i]}) {
		$rrtstdev[$i] += ($val - $rrtmeans[$i]) ** 2;
	}
	$rrtstdev[$i] /= (scalar(@{$rtthash[$i]}) - 1);
	$rrtstdev[$i] = sqrt($rrtstdev[$i]);
}
#The original for loop stupidly caught overrepresented maxishness rather than minishness
#for(my $i=199; $i>-1;$i--) {
for(my $i=3;$i<198;$i++) {
	#Our cutoffs aren't based on smooth data, and so they are vulnerable to local variations in the data.
	#Here we implement a width 3 triangle smoothing function in order to ease the issue
	my $random_window = (($rrtmeans[$i-1] + 7*$rrtstdev[$i-1]) + 2*($rrtmeans[$i] + 7*$rrtstdev[$i]) + ($rrtmeans[$i+1] + 7*$rrtstdev[$i+1]))/4;
	my $sample_window = ($samphash[$i-1] + 2*$samphash[$i] + $samphash[$i+1])/4;
	#Here we implement a width 5 triangle smoothing function in order to ease the issue
#	my $random_window = (($rrtmeans[$i-2] + 7*$rrtstdev[$i-2]) + 2*($rrtmeans[$i-1] + 7*$rrtstdev[$i-1]) + 3*($rrtmeans[$2] + 7*$rrtstdev[$i]) + 2*($rrtmeans[$i+1] + 7*$rrtstdev[$i+1]) + ($rrtmeans[$i+2] + 7*$rrtstdev[$i+2]))/9;
#	my $sample_window = ($samphash[$i-2] + 2*$samphash[$i-1] + 3*$samphash[$i] + 2*$samphash[$i+1] + $samphash[$i+2])/9;
	if($random_window < $sample_window) {
#		print "$i: Mean: " . $rrtmeans[$i] . " + 7*stdev: " . 7*($rrtstdev[$i]) . " = " . $rrtmeans[$i] + 7*($rrtstdev[$i]) . " < " . $samphash[$i] . "\n";
		print "$i: R: $random_window < S: $sample_window\n";
		$cutoff = $i-100;
		print "CUTOFF: " . $cutoff . "\n";
#		print "$cutoff\n";
#		last;
	} elsif ($sample_window > 0){
#		print "LAST: $i: Mean: " . $rrtmeans[$i] . " + 7*stdev: " . 7*$rrtstdev[$i] . " < " . $samphash[$i] . "\n";
		print "$i: R: $random_window < S: $sample_window\n";
		last;
	}
}
